capture log close
log using community.log, replace


*cd "C:\Documents and Settings\agnese\Desktop\communism_1.do"

*========================================================*
*
*  
*
*  First Created: 26/01/09
*
*  Last Edited: 10/02/09 blablabla
*
*  Editor: Agnieszka
*========================================================*

* Data Generation

version 9

drop _all	

set more 1
clear
program drop _all
macro drop _all
matrix drop _all

set virtual on
set mem 500m
set matsize 800
tempfile tmp1 tmp2 tmp3


* Number of observations: check if estimates change greatly if increased....


local obs 10000

set seed 5433   // fix random number generator
set obs `obs'   // number of observations set in local

* generating the U's


gen U0 = invnorm(uniform()) 
gen U1 = 1.2*invnorm(uniform())

* generating the mu's: The ATE = 0 (None in a "random population"

scalar mu1 =  15
scalar mu0 =  15

** Creating Xs

* X1 is the father' years of schooling. 


*gen X1_cont = 8 + 4*invnorm(uniform()) 

*We discretize X1 as follows: 1 = primary school (from 0 to 8); 2 = high school (from 9 to 14) , 3 = university (from 15 onwards)  
 
*gen X1 = X1_cont
*recode X1 min/8 = 0 8/14 = 1 14/max = 2
*tab X1


**X2 is the father belonging to the party

gen X2= (uniform()>=0.8)   
* Draws from uniform and creates 1 if the number is bigger or equal to .8
label var X2 "father belonging to Party"
tab X2



* X3 is S discretized as folllows: 0 if less than 5; 1 if between 5 and 10; 2 if more than 10

gen S = 0.25*(15 + U1) + 0.25*(15 + U0) + 2*invnorm(uniform())
scalar S_bar = 9.358647
gen S_std = S- S_bar

gen X3 = S
recode X3 min/5 = 0 5/10 = 1 10/15 = 2 
label var X3 "social work"
tab X3


* X4 is the daily hours spent studying

*gen X4_cont = 3 + 0.85*invnorm(uniform()) 
* Now we discretize X4 as follows: 0 if less than 1 hour; 1 if between 1 and 2; 2 if between 2 and 3; 3 if between 3 and 4; 4 if more than 4 

*gen X4 = X4_cont
*recode X4 min/1 = 0  1/2 = 1 2/3 = 2 3/4 = 3 4/max = 4

*tab X4
*xi X4
* this drops X4=0

** X5 = 1 is female 
gen X5= round(uniform(),1)
label var X5 "sex"

* now we generate Ys as function of Xs and idiosyncratic gains

gen Y1_m =  3 + 5.8 + 1.53*U1 + 0.092*X2 + 0.08*S + 0.004*X5
*gen Y1_m =  75+20*U1 + 1.2*X2 + X3 + 0.05*X5 - 250

gen Y0_m = 3 + 5.8 + 1.53*U0 + 0.092*X2 + 0.08*S + 0.004*X5

*gen Y0_m =  75+ 20*U0 + 1.2*X2 + X3 + 0.05*X5 - 250

* treatment effects in the population

gen delta = Y1_m - Y0_m


* ==========================================
* different assignment rule: self-selection:
* People do sports if they get net gains...
* eg. they are more competitive...
* ==========================================

gen D_ss = (U1-U0>=0)

* to play around with this, I decided to look at a plot that illustrates our bias a bit... along the lines of ichino's diagrams
* Since in stata i am not sure how to do it easier, here's a really complicated way of doing it :)

gen U1_t = U1 if D_ss ==1
gen U0_t = U0 if D_ss ==1



gen U1_c = U1 if D_ss ==0
gen U0_c = U0 if D_ss ==0

* Comparison 1:*
dotplot U1_t U1_c, title(Idiosyncratic gains from treatment)
pause on
*pause Hit q <RETURN> to continue


dotplot U0_t U0_c, title(Idiosyncratic gains from not doing it)
*pause Hit q <RETURN> to continue


quietly ci U1_t
scalar U1_t_bar = r(mean)


quietly ci U0_t
scalar U0_t_bar = r(mean)


quietly ci U1_c
scalar U1_c_bar = r(mean)

quietly ci U0_c
scalar U0_c_bar = r(mean)



gen TREATED = 1 in 1/2
replace TREATED = U1_t_bar in 1/1
replace TREATED = U0_t_bar in 2/2


gen CONTROL = 1 in 1/2
replace CONTROL = U1_c_bar in 1/2
replace CONTROL = U0_c_bar in 2/2

dotplot TREATED CONTROL, title(Illustration with group means)
*pause Hit q <RETURN> to conclude






***************************************************************** Naive Estimator ************************************************

* Generate the observed sample. Note that we have the counterfactual since we generated them.... 

* Note: multiplying the value Y_ols takes depending on the value of the dummy with the condition on the dummy in brackets!

gen Y_ols = Y0_m*(D_ss==0) + Y1_m*(D_ss==1)

* ==========*
*		*
* Run OLS	*
*		*
* ==========*

reg Y_ols D_ss


* ===================== *
*				*
* IV Estimation 		*
*				*
* ===================== *


gen Z_heli = round(uniform(),1)


* every pupil as his own utility function ...


* selection into treatment with the valid instrument
* a) for a strong instrument (high correlation, i.e. few always- and never-takers

gen D_heli_star = -2+4*Z_heli + (U1 - U0)				

* b) for a weak instrument (low correlation)

gen D_heli_star_weak = -1+0.8*Z_heli + (U1 - U0)


gen D_heli = (D_heli_star >= 0)

gen D_heli_weak = (D_heli_star_weak>=0)

correlate D_heli Z_heli
correlate D_heli_weak Z_heli

* Now find out who is who in this sample.... 

** Alternative Scenario
* Let's create Z_heli_alter as the counterfactual of Z_heli


gen Z_heli_alter = (Z_heli == 0)



gen D_heli_star_alter = -2+4*Z_heli_alter + (U1 - U0)
gen D_heli_alter = (D_heli_star_alter >= 0)

gen D_heli_alter_weak = (-1 + 0.8*Z_heli_alter >= U0-U1)



* To Distinguish the different types of people

* 1. compliers...
* Compliers are those for which the assignment changes the decision to treatment


gen complier = (D_heli == Z_heli & D_heli_alter == Z_heli_alter)  
gen complier_weak = (D_heli_weak == Z_heli & D_heli_alter_weak == Z_heli_alter)  

tab complier 
tab complier_weak


* 2. alwaystakers
* AlwaysTakers are those that always do sport independently of the assignment


gen alwaystaker = (D_heli == 1 & D_heli_alter == 1)
gen alwaystaker_weak = (D_heli_weak == 1 & D_heli_alter_weak == 1)




* 3. nevertakers
* NeverTakers are those that always never do sport independently of the assignment


gen nevertaker = (D_heli == 0 & D_heli_alter == 0 )
gen nevertaker_weak = (D_heli_weak == 0 & D_heli_alter_weak == 0 )

* ================================================================== *
* get an impression about the true ATT for the different types here: *
* ================================================================== *

bysort complier alwaystaker nevertaker: su(delta)
bysort complier_weak alwaystaker_weak nevertaker_weak: su(delta)

bysort complier alwaystaker nevertaker: su(delta) if D_heli == 1
bysort complier_weak alwaystaker_weak nevertaker_weak: su(delta) if D_heli_weak == 1





* Note: we have no alwaystakers in the case of weak instruments ???



gen Y_IV = Y0_m*(D_heli==0) + Y1_m*(D_heli==1)
gen Y_IV_weak = Y0_m*(D_heli_weak==0) + Y1_m*(D_heli_weak==1)


ivregress 2sls Y_IV (D_heli = Z_heli), first
ivregress 2sls Y_IV_weak (D_heli_weak = Z_heli), first

ivreg Y_IV (D_heli = Z_heli), first
ivreg Y_IV_weak (D_heli_weak = Z_heli), first



**IV WITH FULL COMPLIANCE
gen Y_IV_full = Y0_m*(Z_heli==0) + Y1_m*(Z_heli==1)
reg Y_IV_full Z_heli
 
******************************************************
*								     *
*		Regression Discountinuity Design	     *
*    								     *
******************************************************

* S: hours of social work.. 
* those are selected that have more hours of social work because it is appreciated in a communist regime. Note that those who 
* do more work for their country also get better marks at school.

* generate S that is correlated with Y0 and Y1
* regress the Y1, Y0 on S

* Define the cut-off point





gen T = (S>=S_bar)



gen Y_RDD_sharp = Y0_m*(T==0) + Y1_m*(T==1)


* Plot.
twoway (scatter Y_RDD_sharp S), by(T)


rd Y_RDD_sharp S_std, gr bw(1.1)




* (partially) fuzzy RDD *
gen T_fuzzy = T*(delta>=0)
gen Y_RDD_fuzzy = Y0_m*(T_fuzzy==0) + Y1_m*(T_fuzzy==1)


rd Y_RDD_fuzzy T_fuzzy S_std, gr bw(1.1)

* Note: we might have the probability towards high values even go down if we have few observations and those people have idiosyncratic losses
* so they don't select to be treated....

*******************************
*					*
*		MATCHING		*
*******************************

gen V = invnorm(uniform()) // uniform is enough, since we are just interested in value w/ 50% higher / lower than 0.
gen Z_random_ols = 0
replace Z_random_ols = 1 if V>=0


 
*Next step is to create a variable Y_m based on the random assignment to treatment (Z_random_ols). Note here we have full compliance.

*gen Y_m = Y1_m*(Z_random_ols==1)+Y0_m*(Z_random_ols==0)

* Now we match the observations using just X2, X5 and X3

/* We create new variables named box_nw, box_ne; box_se; box_sw as follows
gen box_nw = Y_m if X2 == 0 & X3 == 1 & X5 == 1
gen box_nm = Y_m if X2 == 0 & X3 ==0  & X5 == 1
gen box_ne = Y_m if X2 == 0 & X3 == 1 & X5 == 0
gen box_mw = Y_m if X2 == 1 & X3 == 0 & X5 == 1
gen box_mm = Y_m if X2 == 1 & X3 == 0 & X5 == 0
gen box_me = Y_m if X2 == 1 & X3 == 1 & X5 == 0
gen box_sw = Y_m if X2 == 0 & X3==  0 & X5 == 0
gen box_se = Y_m if X2 == 1 & X5 == 1 & X3 == 1 
gen box_41 = Y_m if X2 == 0 & X3 == 2 & X5 == 0
gen box_42 = Y_m if X2 == 0 & X3 == 2 & X5 == 1
gen box_43 = Y_m if X2 == 1 & X3 == 2 & X5 == 1
gen box_44 = Y_m if X2 == 1 & X3 == 2 & X5 == 0 

* Now we compute the mean for treated and nontreated in each box


su box_ne if Z_random_ols == 1
scalar mean1 = r(mean)
scalar obs1 = r(N)
su box_ne if Z_random_ols == 0
scalar diff1 = mean1 - r(mean)
display "delta box_ne is equal to "diff1

su box_nw if Z_random_ols == 1
scalar mean2 = r(mean)
scalar obs2 = r(N)
su box_nw if Z_random_ols == 0
scalar diff2 = mean2 - r(mean)
display "delta box_nw is equal to "diff2

su box_nm if Z_random_ols == 1
scalar mean3 = r(mean)
scalar obs3 = r(N)
su box_nm if Z_random_ols == 0
scalar diff3 = mean3 - r(mean)
display "delta box_nm is equal to "diff3

su box_mw if Z_random_ols == 1
scalar mean4 = r(mean)
scalar obs4 = r(N)
su box_mw if Z_random_ols == 0
scalar diff4 = mean4 - r(mean)
display "delta box_mw is equal to "diff4

su box_mm if Z_random_ols == 1
scalar mean5 = r(mean)
scalar obs5 = r(N)
su box_mm if Z_random_ols == 0
scalar diff5 = mean5 - r(mean)
display "delta box_mm is equal to "diff5

su box_me if Z_random_ols == 1
scalar mean6 = r(mean)
scalar obs6 = r(N)
su box_me if Z_random_ols == 0
scalar diff6 = mean6 - r(mean)
display "delta box_me is equal to "diff6

su box_se if Z_random_ols == 1
scalar mean7 = r(mean)
scalar obs7 = r(N)
su box_se if Z_random_ols == 0
scalar diff7 = mean7 - r(mean)
display "delta box_se is equal to "diff7

su box_sw if Z_random_ols == 1
scalar mean8 = r(mean)
scalar obs8 = r(N)
su box_sw if Z_random_ols == 0
scalar diff8 = mean8 - r(mean)
display "delta box_sw is equal to "diff8

su box_41 if Z_random_ols ==1
scalar mean9 = r(mean)
scalar obs9 = r(N)
su box_41 if Z_random_ols ==0
scalar diff9 = mean9 - r(mean)
display "delta box_41 is equal to "diff9

su box_42 if Z_random_ols ==1
scalar mean10 = r(mean)
scalar obs10 = r(N)
su box_42 if Z_random_ols ==0
scalar diff10 = mean10 - r(mean)
display "delta box_42 is equal to "diff10

su box_43 if Z_random_ols ==1
scalar mean11 = r(mean)
scalar obs11 = r(N)
su box_43 if Z_random_ols ==0
scalar diff11 = mean11 - r(mean)
display "delta box_43 is equal to "diff11

su box_44 if Z_random_ols ==1
scalar mean12 = r(mean)
scalar obs12 = r(N)
su box_44 if Z_random_ols ==0
scalar diff12 = mean12 - r(mean)
display "delta box_44 is equal to "diff12

*su  box_nw box_nm box_ne box_mw box_mm box_me box_sw box_se box_41 box_42 box_43 box_44


*We take the avarage of the treatment effects in each box and weight them by the number of treated
scalar ATT_random = (diff1*obs1 + diff3*obs3 + diff2*obs2 + diff4*obs4 + diff5*obs5 + diff6*obs6 + diff7*obs7 + diff8*obs8 + diff9*obs9 + diff10*obs10 + diff11*obs11 + diff12*obs12)/(obs1+ obs2+obs3+obs4+obs5+obs6+obs7+obs8+obs9+obs10+obs11+obs12)
display "the ATT_random is equal to "ATT_random */


** Now we want to consider the case of perfect matching but without random assignment to treatment in each box. "Soviet matching"
* to do this we take D_ss and construct Y_ms based on Xs.
gen Y_ms = Y1_m*(D_ss==1)+Y0_m*(D_ss==0)
label var Y_ms "avarage score"

* How serious is the endogeneity problem? correlate Y_ms and D_ss
corr Y_ms D_ss

gen sbox_nw = Y_ms if X2 == 0 & X3 == 1 & X5 == 1
label var sbox_nw "NP1F"
gen sbox_nm = Y_ms if X2 == 0 & X3 ==0  & X5 == 1
label var sbox_nm "NP0F"
gen sbox_ne = Y_ms if X2 == 0 & X3 == 1 & X5 == 0
label var sbox_ne "NP1M"
gen sbox_mw = Y_ms if X2 == 1 & X3 == 0 & X5 == 1
label var sbox_mw "P0F"
gen sbox_mm = Y_ms if X2 == 1 & X3 == 0 & X5 == 0
label var sbox_mm "P0M"
gen sbox_me = Y_ms if X2 == 1 & X3 == 1 & X5 == 0
label var sbox_me "P1M"
gen sbox_sw = Y_ms if X2 == 0 & X3==  0 & X5 == 0
label var sbox_sw "NP0M"
gen sbox_se = Y_ms if X2 == 1 & X5 == 1 & X3 == 1 
label var sbox_se "P1F"
gen sbox_41 = Y_ms if X2 == 0 & X3 == 2 & X5 == 0
label var sbox_41 "NP2M"
gen sbox_42 = Y_ms if X2 == 0 & X3 == 2 & X5 == 1
label var sbox_42 "NP2F"
gen sbox_43 = Y_ms if X2 == 1 & X3 == 2 & X5 == 1
label var sbox_43 "P2F"
gen sbox_44 = Y_ms if X2 == 1 & X3 == 2 & X5 == 0
label var sbox_44 "P2M"

* Now we compute the mean for treated and nontreated within each box

su sbox_ne if D_ss == 1
scalar mean1 = r(mean)
scalar obs1 = r(N)
su sbox_ne if D_ss == 0
scalar diff1 = mean1 - r(mean)
display "delta sbox_ne is equal to "diff1

su sbox_nw if D_ss == 1
scalar mean2 = r(mean)
scalar obs2 = r(N)
su sbox_nw if D_ss == 0
scalar diff2 = mean2 - r(mean)
display "delta sbox_nw is equal to "diff2

su sbox_nm if D_ss == 1
scalar mean3 = r(mean)
scalar obs3 = r(N)
su sbox_nm if D_ss == 0
scalar diff3 = mean3 - r(mean)
display "delta sbox_nm is equal to "diff3

su sbox_mw if D_ss == 1
scalar mean4 = r(mean)
scalar obs4 = r(N)
su sbox_mw if D_ss == 0
scalar diff4 = mean4 - r(mean)
display "delta sbox_mw is equal to "diff4

su sbox_mm if D_ss == 1
scalar mean5 = r(mean)
scalar obs5 = r(N)
su sbox_mm if D_ss == 0
scalar diff5 = mean5 - r(mean)
display "delta sbox_mm is equal to "diff5

su sbox_me if D_ss == 1
scalar mean6 = r(mean)
scalar obs6 = r(N)
su sbox_me if D_ss == 0
scalar diff6 = mean6 - r(mean)
display "delta sbox_me is equal to "diff6

su sbox_se if D_ss == 1
scalar mean7 = r(mean)
scalar obs7 = r(N)
su sbox_se if D_ss == 0
scalar diff7 = mean7 - r(mean)
display "delta sbox_se is equal to "diff7

su sbox_sw if D_ss == 1
scalar mean8 = r(mean)
scalar obs8 = r(N)
su sbox_sw if D_ss == 0
scalar diff8 = mean8 - r(mean)
display "delta sbox_sw is equal to "diff8

su sbox_41 if D_ss == 1
scalar mean9 = r(mean)
scalar obs9 = r(N)
su sbox_41 if D_ss == 0
scalar diff9 = mean9 - r(mean)
display "delta sbox_41 is equal to "diff9

su sbox_42 if D_ss == 1
scalar mean10 = r(mean)
scalar obs10 = r(N)
su sbox_42 if D_ss == 0
scalar diff10 = mean10 - r(mean)
display "delta sbox_42 is equal to "diff10

su sbox_43 if D_ss == 1
scalar mean11 = r(mean)
scalar obs11 = r(N)
su sbox_43 if D_ss == 0
scalar diff11 = mean11 - r(mean)
display "delta sbox_43 is equal to "diff11

su sbox_44 if D_ss == 1
scalar mean12 = r(mean)
scalar obs12 = r(N)
su sbox_44 if D_ss == 0
scalar diff12 = mean12 - r(mean)
display "delta sbox_44 is equal to "diff12

scalar ATT_selection = (diff1*obs1 + diff3*obs3 + diff2*obs2 + diff4*obs4 + diff5*obs5 + diff6*obs6 + diff7*obs7 + diff8*obs8 + diff9*obs9 + diff10*obs10 + diff11*obs11 + diff12*obs12)/(obs1+ obs2+obs3+obs4+obs5+obs6+obs7+obs8+obs9+obs10+obs11+obs12)
display "the ATT_selection is equal to "ATT_selection

*su sbox_nw sbox_nm sbox_ne sbox_mw sbox_mm sbox_me sbox_sw sbox_se sbox_41 sbox_42 sbox_43 sbox_44
tabstat  sbox_nw sbox_nm sbox_ne sbox_mw sbox_mm sbox_me sbox_sw sbox_se sbox_41 sbox_42 sbox_43 sbox_44, by(D_ss) stat(mean, co) nosep


* ========================================= *
* OLS Regression with fully saturated model *
* ========================================= *

* Create Interaction Terms

* For X3 the dummy variables (Reference group: X3 = 0)
tab X3, gen(X3_dum) 

*X2 and X3

gen X2X3_1 = X2*X3_dum2
gen X2X3_2 = X2*X3_dum3

*X2 and X5

gen X2X5 = X2*X5

* X3 and X5

gen X3X5_1 = X5*X3_dum2
gen X3X5_2 = X5*X3_dum3

* all three

gen X2X3X5_1 = X2*X5*X3_dum2
gen X2X3X5_2 = X2*X5*X3_dum3

reg Y_m Z_random_ols X2 X3_dum1 X3_dum2 X5 X2X3_1 X2X3_2 X2X5 X2X3X5_1 X2X3X5_2


reg Y_ms D_ss X2 X3_dum1 X3_dum2 X5 X2X3_1 X2X3_2 X2X5 X2X3X5_1 X2X3X5_2


* ========================= *
* Propensity Score Matching *
* ========================= *


* For the Self-Selection:
quietly pscore D_ss X2 S X5, pscore(p_select) blockid(id_select)

atts Y_ms D_ss, pscore(p_select) blockid(id_select)

*quietly pscore D_ss X2 S X5, pscore(p_select) blockid(id_select) detail

*tab p_select






* This saves the data file....
save  truth.dta, replace



* this closes ALL graphs!!!


*win man close graph _all


log close


view community.log
pause off
